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Abstract 



We develop an efficient algorithm for Synthetic Aperture Sonar imaging based on 
the one-way wave equations. The algorithm utilizes the operator-splitting method 
to integrate the one-way wave equations. The well-posedness of the one-way wave 
equations and the proposed algorithm is shown. A computational result against real 
field data is reported and the resulting image is enhanced by the BV-like regularization. 



1 Introduction 



In this paper we discuss the migration method based on one-way wave equations [H [2] 
for Synthetic Aperture Sonar (SAS) imaging [3j. The one-way wave equation integrates 
the data within a given angle and minimizes the undesirable effects of unwanted reflections. 
Efficient and stable integration methods of the one-way wave equation based on the operator 
splitting method are used to develop a fully discretized algorithm. The stability analysis 



o 
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and the required operation count of the proposed algorithm are given. We test the proposed 
method for real field data and report our SAS imaging results. We also discuss the image 
enhancement method for the resulting images, based on BV-like regularization technique [5]. 

In side-scan (side-looking) sonar systems a platform containing a moderately large real 
aperture antenna travels along a rectilinear path in the along track direction and periodically 
transmits a pulse at an angle that is perpendicular to the platform path. These systems 
produce strip-map (SAS) images . A strip-map image is built up as follows; the imaging 
system operates such that the echoes from the current pulse are received before the next 
pulse is transmitted. As these echoes are received they are demodulated, pulse compressed, 
and detected (only the magnitude information is retained). Each detected pulse produces a 
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range line of the real aperture image. As the platform moves these range lines are displayed 
next to each other at pixel spacings that scale relative to the along track spacing of the 
pulses Ax = v p t where v p is the platform velocity and r is the pulse repetition period. The 
final image is essentially a raster scan of a strip of the sea floor, hence the name " strip-map 
image" . Synthetic aperture imaging is a coherent imaging technique that exploits the extra 
information available in the phase of the real aperture data. We adopt the Stop and Go 
model; a point source radiates at time t = 0, a spherical wave that reaches the sampling 
points after different time intervals. If the source is placed at (xq, Zq) the time t(x, xo, z ) at 
which the wave arrives at the sampling point (x, z) is: 



The field d due to a distribution s(x, z) of source emitting at t = can be expressed by 



along the sonar path (x, z = 0) = T. 

Thus, SAS imaging is formulated as a linear inverse problem; 

Problem: Reconstruct s(x,z) from SAS data SAS(x, t). 

Among a number of algorithms [3j EJ E] and reference therein, which have been developed 
for Problem the frequency domain uj-k method based on Stolt's map [3], [7J [1] is the most 
efficient and accurate method. As will be discussed in Section 4 it has certain limitations, es- 
pecially it assumes the homogeneous scattered media. The proposed method can incorporate 
inhomogeneous media and has additional capabilities, (see Section 4). 

An outline of the paper is as follows. In Section 2 we describe a geometric migration 
method based on one-way wave equations for reconstructing s. A noble algorithm using 
the integration of the one-way wave equations based on the operator-splitting method is 
developed and its stability and complexity are analyzed in Section 3. In Section 3 we list 
advantages of the proposed method comparing to the u-k method. The image enhancement 
technique based on the BV-type reguralization is discussed in Section 4. In Section 5 we 
present a test against real field data, provided by the Naval Surface Warfare Center-Panama 
City, Florida and a comparison with the u-k method. 



t(x,z,x ,z ) = -a/C 



x - x ) 2 + (z - z ) 2 . 




where d is the Fourier transform (in time) of the signal d. SAS measures 



SAS(x,t) = d(x,z = 0,t) 
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2 Geometric Migration 



We construct an approximating solution based on the geometrical migration via the one-way 
wave equations. Let 

A(k x , lo) = T X £ SAS(x, t). 
Assume the plane wave extrapolation 

D(k x , k z , uj) = A(k x , uj)exp(j(ut + k x x + k z z)) 

with 

Then the inverse Fourier transform of D 

1 f 

d(x,z,t) = - — / D(k x , k z ,u)dk x dk z du) 

(27T) J J 

satisfies the wave equation 

4 d 2 d d 2 d d 2 d 



^ c 2 dt 2 dx 2 + dz 2 

with the boundary condition at z = 

d(x,0,t) = SASOM) 

and 

dd 

d(x, z, T) — and z, T) = 0. 

Wave equation based migration integrates the wave equation (3) backward in time to obtain 
an approximation s of distribution s as; 

d(x, z, 0) = s(x, d). 

SAS data is created by integrating over the beam-width of the sensor. The radiation 
pattern of any dimension (width or length) of an aperture has an angular dependence that is 
referred to as the beam pattern of the aperture. Beam patterns are frequency dependent and 

Q 

have beam-widths given by the 3dB response of their main lobes; 6 = a w -j-^ where D is the 

length of the aperture and / are the frequency of the signal that the aperture is transmitting 
or receiving. The term a w is a constant reflecting the main lobe widening due to weighting 
of the aperture illumination function. For example / = 120kHz and D = 0.04m and a w — 1 
gives 9 = 17.9 degrees Thus, in order to speed-up the wave equation based algorithm and 
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minimize the undesirable effects of unwanted reflections we use the (15 degree) one-way wave 
equation based on 



(4) 



where k = 2u/c and we assumed \k x /k\ « 1. In time domain (4) is equivalently written as 
(5) 



4 d 2 u 2 d 2 u 1 d 2 u 
c^W + cdzdi ~ 2dx^' 



with 



u(t,x,0) = SAS(x,t), xeT. 

An advantage of the method is that it allows one to have a specified variable wave speed 
c = c(x, z) of media. The corresponding method for the polar and cylindrical geometry is 
given as 

Polar coordinate 



Cylinder 



4 d 2 u 2 d 2 u 11 d 2 u 
^~dt 2+ c dvdt ~ 2rW' 



4 d 2 u 2 d 2 u 1 ld 2 u d 2 u 



c 2 dt 2 c dvdt 2rd6 2 dz 2 
We can derive the wide angle one-way wave equation based on the rational approximation 



(6) 

we have 



*- = T-(t)^* (1 -i 



a{k x /k) 2 
- (3{k x /kf 



h(k - fa) = k 2 - (a + [3)k 2 x 



The differential form is given by 

4 d 2 u d (2du 



(7) 



c f d 2 u 

r 2 or- ' 777 [~~hJ~ S 2 / JJ^ ( " 



{a + 13) 



d 2 u 
dx 2 



With a = .5, (3 = .25 and a = .478, (3 = .376, (7) is called 45 degree and 65 degree 
approximation, respectively. 



3 Migration by the operator splitting 

Q 

With normalization of the time (t) by the wave speed - and reverting the time, (5) is written 



as 



(8) 




/0 \ 

d 

V° ~TzJ 



u 



/ 1 \ 

i d 2 

V 2&^ ° / 



u 
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So, we apply the time splitting on [t, t + At] of the Lie- Trotter form 



(9) 



/ \ 

d 

V° ~TzJ 



u 



/ 1 \ 

1 d 2 



u 



The first step of (9) is equivalent to the shift operation; 

v(t + At,x,z) = v(t,x,z - At), z>At 
d 

v(t + At, x, z) = — SAS(x, t + z), z E [0, At) 

The second step of (9) is the one-D wave equation in x and is well-posed. In fact, let 

n = \-L,L\ x [0, 1] and H^tQ.) = {(f) E L 2 (fi) : ^-<j> E L 2 (fi)}. Let X 1 = H^tQ.) x L 2 (fi) 

ox 

be the Hilbert space equipped with 

IMI ^ = j/Yx ? + M?)dxdz - 

Define the linear operator A± on X\ by 

1 d 2 u 

Mu,v) = (v,-—) 

with 

dom{A x ) = {(u,v) EX 1 :vE H 1 '*^), E L 2 (tt) 

with tt(±L,z) =0} 
ox 

Then, A\ is dissipative and skew-adjoint on X 1 and thus generates a strongly continuous 
group on X\. Hence, it is easy to show that if (u,v) is generated by (9) then 

rt+At a 

\(u,v)(t + At)\ 2 Xl < \(u,v)(t)\ 2 Xi + J \-SAS(x, s)\ 2 dxds 



and 



\(u,v)(T)\ 2 Xi < / \^SAS(x,t)\ 2 dxdt 



Similarly, we can argue that (8) itself is well-posed, i.e., if we define the operator A on X 1 
by 



A(u,v) = (v,div X:Z (^^-,-v)) 



with 



1 du 
dom(A) = {(u,v) E Xi : v E H 1 '*^), dw x>z (- — , -v ) E L 2 (Q) with v(x,z) = 0, —(±L,z) 

Zi (JJb CJJL 



then A is dissipative and generates a contractive, strongly continuous semigroup on X\. 
We fully discretize (9) and obtain 

Algorithm I 

Offi = v?j, l<j<n with = *— L 

(10) 

u^t 1 = (J + -H)-\u n j + AtuJ 1 ), = — '' J , 1 < j < min(n, M) 

where u" - and i;" - represents the value of u and t> at the grid-point (iAx,jAz) at time nAt, 

Az 

respectively. Here, At = Az and c = — — , H £ r n + 1 > n + 1 j s the tri-diagonal matrix defined 
by 

(ifu)i = -(uj+i - 2«i + tij_i), 2 < i < N, 

and (Hu)i — — {u^ — ui), (Hu)n+i = u^+i — un 

d 2 u 

and corresponds to the central difference approximation of — -7—7. Also, we used the implicit 

dx 2 

Euler scheme to integrate the second step (1-D wave equation in x). That is, 

u? +1 -u? 1 , „ s v n+1 



:(Hu)i, — = u 



n+l 



At Ax 2K /l ' At 

The number of operations at the n-th time step of (10) is of order 0(N min(n, M)). M is 
the number of the focusing step at each pixel in cross-range direction x and if j > M, 
then = u i,j- Thus, the total operation is of order 0(MM 2 ). 

For the wide angle equation (7) we define 



2du c [ 
'cdt~^2j 

It follows from (7) that 



, ' d 2 u , , 2 du 

F=-^-0,l w * and v = - w 



c dt c dt 2 dx 2 dz dx 2 

and 



Thus, (7) is equivalent to 



2 dF dF d 2 u 

+ = Ot ■ 



' 1 1 \ 2 du . „. „ 

,11) -— = (u-F) + F 



c dt dz dx 2 

2du 
c~dt 



2d, „ <9 2 w 
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With v = v — F, we use the three step splitting: 
( dF dF 



'12) 



dt ^ dz 
du 

dv 



du 

m =v 

dv d 2 v 
dt dx 2 



f dF 

~dt 



a 



d 2 F 

dx 2 



du 



dv 
di 







If P — then v — 0, F — v and thus it reduces to the two-step splitting method (9). The 
first equation is accompanied by the boundary condition 

d 

F(t, x, 0) = ^SAS(t, x) - v(t, x, 0). 

Each step of (12) is a well-posed linear system as shown above and we can prove that (11) 
is well-posed. In fact, let Q = [— L, L] x [0, 1] and define the linear operator on A 2 on 
X 2 = L 2 {Q) x H 1 '*^) x L 2 (Q) by 

a ~n / 9F d 2 u „ ^ „d 2 u, 

A(F ,„,„) = ( __ + a _„ + i , /3 _) 



with 



d 

dom(A 2 ) = {-^F e L 2 (tt), with F(-,0) = and 

—u e L 2 (Q) with -(±L, z) = 0, -(v + F) e L 2 m- 



We equip X 2 with norm 



\(F,u,v)\ 2 X2 = [ {\^u\ 2 + -\F\ 2 + -\v\ 2 )dxdz 
Jq ox a p 



Then, A 2 is dissipative, i.e., 

(A 2 (F,u,v),(F,u,v)) 



[,d 2 u,„ du d . _ dF , , 



1 



l)| 2 cfe < 0. 



Since ran^e(^4 2 ) = -^2, -^2 generates a strongly continuous, contraction semigroup on X 2 . 
Similarly, we have the energy estimate 



/ 



(|^ M (T)| 2 + i|F(T)| 2 + iKT)| 2 ^< 



a 



\F(t,x,0)\ 2 dxdt. 



Algorithm I is extended to integrate (12) as follows; 
Algorithm II 

r i,j+l ~ r i,ji 1 — J — n W1XI1 r i,0 — V i,0 

u n + l = (J + f3c 2 H)-\u™ 3 + AtF^ +1 ), F n / 1 = ' J — , 1 <j < min(n, M) 

li"! 1 = (J + ac 2 ^)"^^ 1 + Atw",), < +1 = ' J ' J , \ <j< min(n, M). 
That is, we require double the operations for the integration of Algorithm II. 



4 Advantages of the proposed methods 

The frequency domain ui-k method based on Stolt's map [7] is the most efficient and accurate 
method for the homogeneous media due to the efficiency of fast Fourier transform. It also 
assumes a rectilinear sonar path. 

We can use our proposed algorithms as a means to compensate the motion of sonar path. 
That is, let T be a curved sonar path and T is a reference rectilinear path (z=0). Then we 
solve (5) or (7) on the domain enclosed by the boundaries V and T with boundary value 

u(t, x, z) = SAS(t, x), (x, z) G T 

In this way we have the mapped-SAS data u(t, x, 0) at T and then apply the omega — k 
method for the rectangular domain Q. 

Our implementation (10) of the one-way wave equations is easily adjusted to the case of 
layered media c = c(z) by varying the range increments Az. 

The proposed method can allow to localize the integration on sub-layered regions (assum- 
ing the homogeneous media). Also, we can integrate (5) or (7) in overlapped sub-domains in 
the cross-range (x) direction and then apply the superposition. This improves the efficiency 
of the proposed algorithms. 



5 BV-type Regular izat ion for Enhancement of S AS imag- 
ing 

SAS imaging s(x, z) may be altered by inhomogeneity of the field, sensor noise and irreg- 
ularity of the sonar path and so on. We use the image enhancement technique based on 
BV-type reguralization [5]. 



8 



Enhancement S of s minimizes 

I 2 d-rdv. 4- R I m(\. 

dx dz 
(3 > is the regularization parameter 

<M = / *l VSR deBne s the restoration energy. 



(12) / \S - s\ 2 dxdz + (3 / <^(|^| 2 + |^| 2 )cte^ 

where 



and 



The followings summarize our findings in |5] on the enhancement based on(12); 

• y?(£ 2 ) = t 2 corresponds to the standard Gaussian filter and works well for a smooth 
image s. 

• ip{t 2 ) = t corresponds to the BV (nonlinear) filter and restores edges and flat regions 
of image s very well. But, it has significant stair-case effects. 

• In order to deal with images with multi-scales of edges, flat, and smooth regions we 
developed an algorithm which uses 

( 1 



_ s e l,oo 
1 s G [5, 1] 
4= ^(0,5) 



It is based on the scale analysis and we demonstrated the applicability of the algorithm 
in [5]. 

• The necessary and sufficient condition of (12) is given by 

-/3V- ( V \\VS\ 2 )VS) + S = s. 

An efficient algorithm for finding S based on the fixed point iterate; 

-P V • ( V '(\VS k \ 2 )VS k+1 ) + S k+1 = s 
is developed and analyzed in [5] and is used in our test. 

6 A Test 

The algorithm is successfully applied to real data that are available to us via the Naval Surface 
Warfare Center (NSWC) and shows a promising capability. A full capability is going to be 
tested in the line of its advantages discussed in Section 4. In a CRSC tereport, CRSC- 
TR09-12 at http: / / www.ncsu.edu / crsc / reports / reports09.htm| we show the raw SAS data, 
SAS imaging by algorithm (10), and the image enhanced by our enhancement algorithm. 
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